/*	
	This file is part of the Snoopy code.

    Snoopy code is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    Snoopy code is distributed in the hope that it will be useful,
    but WITHOUT Ag.nyrtot WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Snoopy code.  If not, see <http://www.gnu.org/licenses/>.
*/



#include "snoopy.h"

#include <stdlib.h>

// Timing support
#ifndef MPI_SUPPORT
#ifndef _OPENMP
#include <stdio.h>
#include <time.h>
#endif
#endif

#include "error.h"
#include "gfft.h"
#include "debug.h"

/* shared grid over the entire code */
struct Grid g;

double	*wr1,	*wr2,	*wr3;		/** Temporary real array (alias of complex w**) */
double	*wr4,	*wr5,	*wr6;		/** Temporary real array (alias of complex w**) */
double	*wr7,	*wr8,	*wr9;		/** Temporary real array (alias of complex w**) */
double	*wr10,	*wr11,	*wr12;
double	*wr13,	*wr14,	*wr15;

double complex		*w1,	*w2,	*w3;	/**< Temporary complex array (alias of real wr**) */
double complex		*w4,	*w5,	*w6;	/**< Temporary complex array (alias of real wr**) */
double complex		*w7,	*w8,	*w9;	/**< Temporary complex array (alias of real wr**) */
double complex		*w10,	*w11,	*w12;
double complex		*w13,	*w14,	*w15;



// Parameters
struct Parameters			param;

// Physics variables 
double	nu;
#ifdef BOUSSINESQ
double	nu_th;
#endif
#ifdef MHD
double	eta;
#endif

int		nthreads;								/**< Number of OpenMP threads */


/* Function prototypes */
void init_grid(int nx, int ny, int nz, struct Grid *gridin);
void allocate_field(struct Field *fldi);
void deallocate_field(struct Field *fldi);
void init_N2_profile();
void init_real_mask();

/***************************************************************/
/**
	Init all global variables, aligning them in memory
*/
/***************************************************************/
void init_common(void) {
	/* This routine will initialize everything */
	int i,j,k;
	
	DEBUG_START_FUNC;
	
#ifdef MPI_SUPPORT	
	fftw_mpi_init();
#endif
#ifdef _OPENMP
	if( !(fftw_init_threads()) ) ERROR_HANDLER( ERROR_CRITICAL, "Threads initialisation failed");
#endif
	
    // Init the grid g
    init_grid(param.nx, param.ny, param.nz, &g);
    
    
// Allocate fields
// Complex fields

	w1 = (double complex *) fftw_malloc( sizeof(double complex) * g.nloc);
	if (w1 == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for w1 allocation");
	
	w2 = (double complex *) fftw_malloc( sizeof(double complex) * g.nloc);
	if (w2 == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for w2 allocation");
	
	w3 = (double complex *) fftw_malloc( sizeof(double complex) * g.nloc);
	if (w3 == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for w3 allocation");
	
	w4 = (double complex *) fftw_malloc( sizeof(double complex) * g.nloc);
	if (w4 == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for w4 allocation");
	
	w5 = (double complex *) fftw_malloc( sizeof(double complex) * g.nloc);
	if (w5 == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for w5 allocation");
	
	w6 = (double complex *) fftw_malloc( sizeof(double complex) * g.nloc);
	if (w6 == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for w6 allocation");
	
	w7 = (double complex *) fftw_malloc( sizeof(double complex) * g.nloc);
	if (w7 == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for w7 allocation");
	
	w8 = (double complex *) fftw_malloc( sizeof(double complex) * g.nloc);
	if (w8 == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for w8 allocation");
	
	w9 = (double complex *) fftw_malloc( sizeof(double complex) * g.nloc);
	if (w9 == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for w9 allocation");
	
	w10 = (double complex *) fftw_malloc( sizeof(double complex) * g.nloc);
	if (w10 == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for w10 allocation");
	
	w11 = (double complex *) fftw_malloc( sizeof(double complex) * g.nloc);
	if (w11 == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for w11 allocation");
	
	w12 = (double complex *) fftw_malloc( sizeof(double complex) * g.nloc);
	if (w12 == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for w12 allocation");
	
	w13 = (double complex *) fftw_malloc( sizeof(double complex) * g.nloc);
	if (w13 == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for w13 allocation");
	
	w14 = (double complex *) fftw_malloc( sizeof(double complex) * g.nloc);
	if (w14 == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for w14 allocation");
	
	w15 = (double complex *) fftw_malloc( sizeof(double complex) * g.nloc);
	if (w15 == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for w15 allocation");
			
	/* Will use the same memory space for real and complex fields */
	
	wr1 = (double *) w1;
	wr2 = (double *) w2;
	wr3 = (double *) w3;
	wr4 = (double *) w4;
	wr5 = (double *) w5;
	wr6 = (double *) w6;
	wr7 = (double *) w7;
	wr8 = (double *) w8;
	wr9 = (double *) w9;
	wr10 = (double *) w10;
	wr11 = (double *) w11;
	wr12 = (double *) w12;
	wr13 = (double *) w13;
	wr14 = (double *) w14;
	wr15 = (double *) w15;

	
// Physic initialisation
//	init_real_mask();
	
	nu = 1.0 / param.reynolds;
#ifdef BOUSSINESQ	
	nu_th = 1.0 / param.reynolds_th;
#endif
#ifdef MHD
	eta = 1.0 / param.reynolds_m;
#endif
	DEBUG_END_FUNC;
	return;
}

void finish_common(void) {
	
	fftw_free(w1);
	fftw_free(w2);
	fftw_free(w3);
	fftw_free(w4);
	fftw_free(w5);
	fftw_free(w6);
	fftw_free(w7);
	fftw_free(w8);
	fftw_free(w9);
	fftw_free(w10);
	fftw_free(w11);
	fftw_free(w12);
	fftw_free(w13);
	fftw_free(w14);
	fftw_free(w15);

	return;
}

/*********************************************/
/**
 Create a grid structure with nx, ny nz points
 and decompose over Nprocx, Nprocy procs
 
 @param nx: number of points in the x direction
 @param ny: number of points in the y direction
 @param nz: number of points in the z direction
 @param Nprocx: number of cores in the x direction
 @param Nprocy: number of cores in the y direction
 **/
/*********************************************/

void init_grid(int nx, int ny, int nz, struct Grid *gridin) {
    int nproc;
    int i,j,k;
    ptrdiff_t temp1,temp2;
    DEBUG_START_FUNC;
    
    /* Initialize the grid structure */
    gridin->nxrtot=nx;
    gridin->nyrtot=ny;
    gridin->nzrtot=nz;
#ifdef WITH_2D
	gridin->nzrtot=1;
#endif
    
    gridin->nrtot = nx*ny*nz;
    
    /* number of Complex points */
#ifndef WITH_2D
    gridin->nxtot=gridin->nxrtot;
    gridin->nytot=gridin->nyrtot;
    gridin->nztot=gridin->nzrtot/2+1;
#else
	gridin->nxtot=gridin->nxrtot;
	gridin->nytot=gridin->nyrtot/2+1;
	gridin->nztot=1;
#endif
    
    gridin->ntot = gridin->nxtot*gridin->nytot*gridin->nztot;
    
    /* number of running cores */
#ifdef MPI_SUPPORT
    MPI_Comm_size(MPI_COMM_WORLD, &nproc);
#else
    nproc = 1;
#endif
    
    gridin->nProcx=nproc;
    gridin->nProcy=1;
    
    gridin->nProc = gridin->nProcx*gridin->nProcy;
    
    /*Domain decomposition*/
#ifdef MPI_SUPPORT
    MPI_Comm_rank(MPI_COMM_WORLD,&gridin->rankx);
    if(gridin->rankx==0) gridin->main=1;
    else gridin->main=0;
#else
    gridin->rankx=0;
    gridin->main=1;
#endif
    gridin->ranky=0;
    
#ifdef MPI_SUPPORT
    gridin->nloc = fftw_mpi_local_size_3d(gridin->nxtot, gridin->nytot, gridin->nztot,MPI_COMM_WORLD, &temp1, &temp2);
    /* These are required since local_size_3D returns ptrdiff_t which are machine dependent */
    gridin->nxloc=(int)temp1;
    gridin->startx=(int)temp2;
#else
    gridin->nloc=gridin->ntot;
    gridin->nxloc=gridin->nxtot;
    gridin->startx=0;
#endif
    gridin->nyloc=gridin->nytot;
    gridin->nzloc=gridin->nztot;
    gridin->starty = 0;
    gridin->startz = 0;
    
    /* Map real number array on the complex array */
    gridin->nrloc=2*gridin->nloc;
	gridin->nxrloc=gridin->nxloc;
	gridin->nyrloc=gridin->nyloc;
	gridin->nzrloc=2*gridin->nzloc;
	
    MPI_Printf("init_grid: (nxr, nyr, nzr)   =   (%d,%d,%d)\n",gridin->nxrtot,gridin->nyrtot,gridin->nzrtot);
    MPI_Printf("init_grid: (nProcx, nProcy)  =   (%d,%d)\n",gridin->nProcx,gridin->nProcy);
    MPI_Printf("init_grid: (nxloc, nyloc, nzloc)=(%d, %d, %d)\n",gridin->nxloc,gridin->nyloc,gridin->nzloc);
    
    /* Initialise Fourier coordinate system */
	gridin->kx = (double *) fftw_malloc( sizeof(double) * gridin->nloc );
	if (gridin->kx == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for kx allocation");
	
	gridin->ky = (double *) fftw_malloc( sizeof(double) * gridin->nloc );
	if (gridin->ky == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for ky allocation");
	
	gridin->kz = (double *) fftw_malloc( sizeof(double) * gridin->nloc );
	if (gridin->kz == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for kz allocation");
	
	gridin->kxt = (double *) fftw_malloc( sizeof(double) * gridin->nloc );
	if (gridin->kxt == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for kxt allocation");
	
	gridin->k2t = (double *) fftw_malloc( sizeof(double) * gridin->nloc );
	if (gridin->k2t == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for k2t allocation");
	
	gridin->ik2t = (double *) fftw_malloc( sizeof(double) * gridin->nloc );
	if (gridin->ik2t == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for ik2t allocation");
    
    
	for( i = 0; i < gridin->nxloc; i++) {
		for( j = 0; j < gridin->nyloc; j++) {
			for( k = 0; k < gridin->nzloc; k++) {
				gridin->kx[ IDX3D ] = (2.0 * M_PI) / param.lx * ( fmod(gridin->startx + i + gridin->nxtot /2, gridin->nxtot) - gridin->nxtot / 2);
                
#ifdef WITH_2D
                gridin->ky[ IDX3D ] = (2.0 * M_PI) / param.ly * j;
                gridin->kz[ IDX3D ] = 0.0;
#else
                gridin->ky[ IDX3D ] = (2.0 * M_PI) / param.ly * ( fmod(gridin->starty + j + gridin->nytot /2, gridin->nytot) - gridin->nytot / 2);
				gridin->kz[ IDX3D ] = (2.0 * M_PI) / param.lz * k;
#endif
                
				gridin->kxt[ IDX3D ]= gridin->kx[IDX3D];
                
				gridin->k2t[ IDX3D ] = gridin->kxt[IDX3D] * gridin->kxt[IDX3D] +
                                       gridin->ky[IDX3D]  * gridin->ky[IDX3D] +
                                       gridin->kz[IDX3D]  * gridin->kz[IDX3D];
                
				if ( gridin->k2t[IDX3D] == 0.0 ) gridin->ik2t[IDX3D] = 1.0;
				else	gridin->ik2t[IDX3D] = 1.0 / gridin->k2t[IDX3D];
			}
		}
	}
	
	gridin->kxmax = 2.0 * M_PI/ param.lx * ( (gridin->nxrtot / 2) - 1);
	gridin->kymax = 2.0 * M_PI/ param.ly * ( (gridin->nyrtot / 2) - 1);
	gridin->kzmax = 2.0 * M_PI/ param.lz * ( (gridin->nzrtot / 2) - 1);
#ifdef WITH_2D
	gridin->kzmax = 0.0;
#endif
	
	gridin->kmax=pow(gridin->kxmax*gridin->kxmax+gridin->kymax*gridin->kymax+gridin->kzmax*gridin->kzmax,0.5);
	
	/* Initialize the dealiazing mask Or the nyquist frequency mask (in case dealiasing is not required) */
	
	gridin->mask = (double *) fftw_malloc( sizeof(double) * gridin->nloc );
	if (gridin->mask == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for mask allocation");
	
    for( i = 0; i < gridin->nxloc; i++) {
		for( j = 0; j < gridin->nyloc; j++) {
			for( k = 0; k < gridin->nzloc; k++) {
                
				gridin->mask[ IDX3D ] = 1.0;
				if(param.antialiasing) {
					if( fabs( gridin->kx[ IDX3D ] ) > 2.0/3.0 * gridin->kxmax)
						gridin->mask[ IDX3D ] = 0.0;
                    
					if( fabs( gridin->ky[ IDX3D ] ) > 2.0/3.0 * gridin->kymax)
						gridin->mask[ IDX3D ] = 0.0;
#ifndef WITH_2D
					if( fabs( gridin->kz[ IDX3D ] ) > 2.0/3.0 * gridin->kzmax)
						gridin->mask[ IDX3D ] = 0.0;
#endif
				}
				else {
					if (  gridin->startx + i == gridin->nxtot / 2 )
						gridin->mask[ IDX3D ] = 0.0;
					if ( j == gridin->nytot / 2 )
						gridin->mask[ IDX3D ] = 0.0;
#ifndef WITH_2D
					if ( k == gridin->nztot )
						gridin->mask[ IDX3D ] = 0.0;
#endif
				}
			}
		}
	}
    
	if(param.antialiasing) {
		gridin->kxmax = gridin->kxmax * 2.0 / 3.0;
		gridin->kymax = gridin->kymax * 2.0 / 3.0;
		gridin->kzmax = gridin->kzmax * 2.0 / 3.0;
		gridin->kmax = gridin->kmax * 2.0 / 3.0;
	}
    
}
/*********************************************/
/**
	Allocate a field structure according to the code
	current configuration
	This routine allows one to add extra fields
	to the code very easily.

	@param *fldi: pointer to a field structure to initialize
**/
/*********************************************/
void allocate_field(struct Field *fldi) {
	int current_field, i;
	
	DEBUG_START_FUNC;
	
	// We want to allocate a field structure
	fldi->nfield = 3;	
		
#ifdef BOUSSINESQ
	fldi->nfield++;
#endif

#ifdef MHD
	fldi->nfield=fldi->nfield+3;
#endif

#ifdef COMPRESSIBLE
	fldi->nfield=fldi->nfield+1;
#endif

	// Now we want to initialize the pointers of the field structure
	
	// farray will point to each of the array previously allocated
	fldi->farray = (double complex **) fftw_malloc( sizeof(double complex *) * fldi->nfield);
	if (fldi->farray == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for fldi->farray allocation");
	
	// fname will point to the name of each field
	fldi->fname = (char **) fftw_malloc(sizeof(char *) * fldi->nfield);
	if (fldi->fname == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for fldi->fname allocation");
	
	// Initialise the pointers
	for(i=0 ; i < fldi->nfield ; i++) {
		fldi->fname[i] = (char *) fftw_malloc(sizeof(char) * 10); // 10 character to describe each field
		if (fldi->fname[i] == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for fldi->fname[i] allocation");
	}
	
	// Allocate the arrays and put the right value in each pointer
	
	current_field = 0;

	fldi->vx = (double complex *) fftw_malloc( sizeof(double complex) * g.nloc);
	if (fldi->vx == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for fldi->vx allocation");
	fldi->farray[current_field] = fldi->vx;
#ifndef COMPRESSIBLE
	sprintf(fldi->fname[current_field],"vx");
#else
	sprintf(fldi->fname[current_field],"px");
#endif
	current_field++;
	
	fldi->vy = (double complex *) fftw_malloc( sizeof(double complex) * g.nloc);
	if (fldi->vy == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for fldi->vy allocation");
	fldi->farray[current_field] = fldi->vy;
#ifndef COMPRESSIBLE
	sprintf(fldi->fname[current_field],"vy");
#else
	sprintf(fldi->fname[current_field],"py");
#endif
	current_field++;
	
	fldi->vz = (double complex *) fftw_malloc( sizeof(double complex) * g.nloc);
	if (fldi->vz == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for fldi->vz allocation");
	fldi->farray[current_field] = fldi->vz;
#ifndef COMPRESSIBLE
	sprintf(fldi->fname[current_field],"vz");
#else
	sprintf(fldi->fname[current_field],"pz");
#endif
	current_field++;
	
#ifdef BOUSSINESQ
	fldi->th = (double complex *) fftw_malloc( sizeof(double complex) * g.nloc);
	if (fldi->th == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for fldi->th allocation");
	fldi->farray[current_field] = fldi->th;
	sprintf(fldi->fname[current_field],"th");
	current_field++;
#endif
#ifdef MHD
	fldi->bx = (double complex *) fftw_malloc( sizeof(double complex) * g.nloc);
	if (fldi->bx == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for fldi->bx allocation");
	fldi->farray[current_field] = fldi->bx;
	sprintf(fldi->fname[current_field],"bx");
	current_field++;
	
	fldi->by = (double complex *) fftw_malloc( sizeof(double complex) * g.nloc);
	if (fldi->by == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for fldi->by allocation");
	fldi->farray[current_field] = fldi->by;
	sprintf(fldi->fname[current_field],"by");
	current_field++;
	
	fldi->bz = (double complex *) fftw_malloc( sizeof(double complex) * g.nloc);
	if (fldi->bz == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for fldi->bz allocation");
	fldi->farray[current_field] = fldi->bz;
	sprintf(fldi->fname[current_field],"bz");
	current_field++;
#endif
#ifdef COMPRESSIBLE
	fldi->d = (double complex *) fftw_malloc( sizeof(double complex) * g.nloc);
	if (fldi->d == NULL) ERROR_HANDLER( ERROR_CRITICAL, "No memory for fldi->d allocation");
	fldi->farray[current_field] = fldi->d;
	sprintf(fldi->fname[current_field]," d");
	current_field++;
#endif
	// Add a field here if you need one... (don't forget to ajust fldi.nfield accordingly)
	// *


#ifdef WITH_PARTICLES	
	// Init space for particle storage
	fldi->part = (struct Particle *) malloc(sizeof(struct Particle) * param.particles_n);
	
#endif
	// Ok, all done...
	
	DEBUG_END_FUNC;
	
	return;
}
/*********************************************/
/**
Deallocate a field structure created by
allocate_field
**/
/*********************************************/
void deallocate_field(struct Field *fldi) {
	int i;
	// Free a field structure
	
	DEBUG_START_FUNC;
	
	for(i=0 ; i < fldi->nfield ; i++) {
		fftw_free(fldi->fname[i]);
		fftw_free(fldi->farray[i]);
	}
	
	fftw_free(fldi->farray);
	fftw_free(fldi->fname);
	
#ifdef WITH_PARTICLES
	free(fldi->part);
#endif
	
	// Done
	
	DEBUG_END_FUNC;
	
	return;
}
	

/*********************************************/
/**
Customized random number generator
Allow one to have consistant random numbers
generators on different architectures.
**/
/*********************************************/
double randm(void) {
	const int a	=	16807;
	const int m =	2147483647;
	static int in0 = 13763;
	int q;
	
	// When using mpi, this allows us to have different number series in each process...
	if(in0 == 13763) in0 += 2543 * (g.ranky+g.nProcy*g.rankx);
	
	/* find random number  */
	q= (int) fmod((double) a * in0, m);
	in0=q;
	
	return((double)q/(double)m);
}
/*********************************************/
/**
	 * Normal distribution
	 * Algorithm by D.E. Knut, 1997, The Art of Computer Programmin, Addison-Wesley. 
	 */
/*********************************************/
	 
double randm_normal(void) {
	double v1, v2;
	double rsq=1.0;
	
	while(rsq>=1. || rsq==0.0) {
		v1=2.*randm()-1.0;
		v2=2.*randm()-1.0;
		rsq=v1*v1+v2*v2;
	}
	
	return( v1*sqrt(-2.0 * log(rsq) / rsq));
}

/****************************************************/
/**
	Remove the divergence from a 3D field using
	the projector operator:
	q=q-k.q/k^2
	
	@param qx: x component of the field
	@param qy: y component of the field
	@param qz: z component of the field
*/
/****************************************************/
	
void projector( double complex qx[],
			    double complex qy[],
			    double complex qz[]) {
				
	int i;
	double complex q0;
	
	DEBUG_START_FUNC;
	
	for( i = 0 ; i < g.nloc ; i++) {
		q0 = g.kxt[i] * qx[i] + g.ky[i] * qy[i] + g.kz[i] * qz[i];
		qx[i] = qx[i] - g.kxt[i] * q0 * g.ik2t[i];
		qy[i] = qy[i] - g.ky[i]  * q0 * g.ik2t[i];
		qz[i] = qz[i] - g.kz[i]  * q0 * g.ik2t[i];
	}
	
	DEBUG_END_FUNC;
	
	return;
}

/*********************************************/
/** Compute the energy of a given field.
	@param q complex array containing the field for which
				  we want the total energy
*/
/*********************************************/

double energy(const double complex q[]) {
	
	int i,j,k;
	double energ_tot;
	
	energ_tot=0.0;
	
	for( i = 0; i < g.nxloc; i++) {
		for( j = 0; j < g.nyloc; j++) {
			for( k=0; k < g.nzloc; k++) {
#ifdef WITH_2D
				if( j == 0)
#else
				if( k == 0)
#endif
					// k=0, we have all the modes.
					energ_tot = energ_tot + creal( 0.5 * q[ IDX3D ] * conj( q[ IDX3D ] ) ) / ((double) g.nrtot*g.nrtot);
				else
					// k>0, only half of the complex plane is represented.
					energ_tot = energ_tot + creal( q[ IDX3D ] * conj( q[ IDX3D ] ) ) / ((double) g.nrtot*g.nrtot);
			}
		}
	}
//	energ_tot = 0;
	return(energ_tot);
}

/********************************************/
/**
Return the localtime in seconds. Use different
implementation depending on the avaiable
libraries  **/
/********************************************/
double get_c_time(void) {
#ifdef MPI_SUPPORT
	// We have MPI
	return(MPI_Wtime());
#else
#ifdef _OPENMP
	// We don't have MPI, but we have OpenMP
	return(omp_get_wtime());
#else
	// We really have nothing...
	clock_t now;
	now = clock();
	
	return( (double) now / ( (double) CLOCKS_PER_SEC));

#endif
#endif
}

/******************************************/
/**
	Reduce a variable over all the avaiable processes
	Can add a value on all the process, find a maximum
	or a minimum.
	
	NOTE: This routine makes sense only when MPI_SUPPORT
	is set. If not, this routine does nothing.
	
	@param *var: variable to be reduced
	@param op: operation needed to be done. Can be set to:
		1= Sum over all the processes
		2= Find the maximum over all the processes
		3= Find the minimum over all the processes
*/
/*******************************************/
	
void reduce(double *var, const int op) {
	// op=1 ADD
	// op=2 Max
	// op=3 Min
	
#ifdef MPI_SUPPORT
	double mpi_temp;
	
	mpi_temp=*var;
	
	if(op==1) MPI_Allreduce( &mpi_temp, var, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
	if(op==2) MPI_Allreduce( &mpi_temp, var, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
	if(op==3) MPI_Allreduce( &mpi_temp, var, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
	
#endif	

	// If no MPI, then this routine does nothing...
	return;
}

/* ****************************************************************************/
/** Determines if the machine is little-endian.  If so, 
    it will force the data to be big-endian. 
	@param in_number floating point number to be converted in big endian */
/* *************************************************************************** */

float big_endian(float in_number)
{
    static int doneTest = 0;
    static int shouldSwap = 0;
	
    if (!doneTest)
    {
        int tmp1 = 1;
        unsigned char *tmp2 = (unsigned char *) &tmp1;
        if (*tmp2 != 0)
            shouldSwap = 1;
        doneTest = 1;
    }

    if (shouldSwap)
    {
		unsigned char *bytes = (unsigned char*) &in_number;
        unsigned char tmp = bytes[0];
        bytes[0] = bytes[3];
        bytes[3] = tmp;
        tmp = bytes[1];
        bytes[1] = bytes[2];
        bytes[2] = tmp;
    }
	return(in_number);
}

/******************************************************************************/
/** Test a double for Not a Number error 
	@param xi double to be checked */
/* ************************************************************************** */

void c_nan(double xi, const char ErrorRoutine[], const int line, const char Filename[]) {
	if(isnan(xi)) {
		error_h( 3, "Not a number detected", ErrorRoutine, line, Filename);
	}
	return;
}

#ifdef COMPRESSIBLE
/***************************************************************/
/** Check that the field is definite positive *******************/
/****************************************************************/

void check_positivity(double *wri) {
	int i;
#ifdef _OPENMP
	#pragma omp parallel for private(i) schedule(static)	
#endif
	for(i=0 ; i < g.nrloc ; i++) {
		if(wri[i] < 1.0e-2*g.nrtot) wri[i] = 1.0e-2*g.nrtot;
	}
}
#endif
